close all
clear all
set(0,'DefaultTextFontSize',25)
set(0,'DefaultAxesFontSize',25)
set(0,'DefaultLineLineWidth',4)
set(0,'DefaultAxesLineWidth',4)
set(0,'DefaultSurfaceLineWidth',2)
disp('^^^^^^^^^^^^^^^^^')
tic

%MS notes
% It has now come to be recognized that themanagement of internationally
% shared fishery resources is one of the most significant, if not the most significant,
% resource management issues to have arisen under the EEZ regime.
% The FAO estimates that as much as one-third of world marine capture
% fishery harvests are based upon such shared stocks (Munro et al., 2004).(Munro 2008)


% 
% 
% global theta  numpatches Dii   domain_suit_hab max_age age_mature...
%     age_legal  Wij Rmax alpha1 delta Kvirgin  pr max_simulation_length
%     ...
%        Fsum  hide0_show1_figures K Bi_legal_fish  beta_i_initial ...
%        runTURForder TURFprofitshare Bij_initial TargetBiomass CRgoal theta_percent...

global theta  numpatches Dii   domain_suit_hab max_age age_mature...
    age_legal  Wij Rmax alpha1 delta Kvirgin  pr max_simulation_length ...
       Fsum  hide0_show1_figures K Bi_legal_fish  beta_i_initial ...
       Bij_initial TargetBiomass CRgoal theta_percent...
       FleetProfit_Fsum_opt MPA_locs Payoff_Conservation_EF_alpha_dynamic...
       x0_Fsumguess  LB UB Fsum_UB_prop_wrt_sumFishingHab...
       Bi_initial Bi_Fsum_opt Profit_i_Fsumopt_NoMPAs Payoff_Conservation_EF_alpha_STATIC1...
       NCI_MPA_locs Kunfished_hat...
       NPVon discounted_value_iy discounted_value_ijy
   
%CHOOOOOOOOOOOOOOOOOOOOSE
% max_simulation_length=300;
max_simulation_length=100; %SHORTENED!
plot_parameters=1; 
tune_parameters_1YES_0NO=1;
Tune_iters=3; %5
Fsum_NumIters=10; %pub quality=100
stop_before_optimization=0; %1=yes
Fsum_UB_prop_wrt_sumFishingHab=0.1; %see Fsum_UB_wrt_fishinghabitat.xlsx
prop_MPAs=[0:0.5:1]; %always include 0 and 1. Note yet used
% path('\\tsclient\H\Data_Models\npsol;',path); %This enables npsol
% path('H:\Data_Models\MATLAB\npsol;',path);
% path('C:\Documents and Settings\Dan\Desktop\Crow_SFG1\npsol;',path);
% path('C:\Users\crow\Desktop\npsol;',path);
% path('C:\Documents and Settings\crow\Desktop\Crow\npsol;',path); %SFG2
% path('C:\Documents and Settings\crow\Desktop\Crow SFG3\npsol;',path);
% path('C:\Users\cwhite31\Documents\Data_Models\MATLAB\npsol;',path);
reset_rand_0no_1yes=0;
if reset_rand_0no_1yes==1
    rand('state',sum(1e5.*clock));
    randn('state',sum(1e5.*clock));
end
NPVon=0; %LEAVE (changes later when go to Net present value)
    

% load spatial data
existMPAs=unique([63 68 73 77 79 82 88 89 94 96 135]'); %NCI existing MPAs
load urchinHab135; %area rocky habitat 0-100 m depth
domain_suit_hab=urchinHab135;
MPA_locs=zeros(1,length(domain_suit_hab)); %initially no MPAs
load connect_RedSeaUrchin;
Dii=connect;
% Dii=rand(size(connect)); %random
% Dii=diag(ones(length(connect(:,1)),1)); %local dispersal only
%  Dii=ones(size(connect)); %larval pool
%row-sum standardize
% Diirowsum = sum(Dii,2);
% for tmp=1:length(Dii(:,1)); Dii(tmp,:)=Dii(tmp,:)/Diirowsum(tmp); end
% Diirowsum = sum(Dii,2); %should be one
% if unique(Diirowsum) ~= 1
%     disp('Dii not standardized. Break')
%     break
% end
% for d=1:length(Dii(1,:)); Dii(d,d)=0; end %set diag to zero NO SELF-RECRUITMENT

load Patches
figure
plot(site.Lon,site.Lat, 'bo'); hold on
for tmp=1:length(site.Lon); text(site.Lon(tmp),site.Lat(tmp),num2str(tmp)); end
TargetBiomass = 1e10; %kg 1e6????
numpatches=size(domain_suit_hab,1);   
%ages
age_mature=4; %the age of first reproduction in years
age_legal=8; %the age at which individuals are avaliable to be fished
max_age=20; %maximum age the species reaches
delta=0.08; %The instantaneous rate of non-fishing mortality, expressed in per year units
%von Bertalanffy parameters (L in cm TL, t in years)
Loo=11;
k=0.22;
t0=0;
Lj=Loo.*(1-exp(-k.*((0:max_age)-t0))); %length at each age class
%Weight-at-length equation, W = C1*L^(C2), L in cmTL, W in kg
C1=0.000676304;
C2=2.68;
Wj=C1.*(Lj.^C2); %weight at each age class
%density dependence
%A measure of the strength of density dependence.  
%The ratio of maximum settler survival rate to settler survival rate in the
%unfished system
CRgoal=4;
%guesses
% Rmax = 1412.2422; alpha1 = 1.8723; %targetbiomass=1e6
% Rmax = 141224.2207; alpha1 = 1.8723; %targetbiomass=1e8
Rmax = 1486.8695; alpha1 = 54.8165; %Dii not standardize; targetbiomass=1e6
Rmax = 148686.95; alpha1 = 54.8165; %Dii not standardize; targetbiomass=1e8
Rmax = 14868695.398; alpha1 = 54.8165; %Dii not standardize; targetbiomass=1e10
warning off MATLAB:divideByZero
beta_i_initial=alpha1./(Rmax.*domain_suit_hab); %related to max biomass of recruits in each cell
warning on MATLAB:divideByZero
% ECONOMICS
%MLPA Initiative Draft Regional Profile from 2008 (attached), 
% which indicates that average annual $/lb for the Santa Barbara County 
% sea urchin fishery (1998-2007) is $0.80/lb. ($3,451,281/4,299,401 lb.). 
kg_per_lb=0.45359237;
pr=0.8*kg_per_lb; %$/kg
theta_percent=0.1; %percent of Carrying capacity below which marginal cost=marginal revenue

if plot_parameters==1
figure; surface(Dii); title('Larval dispersal'); xlabel('Destination'); ylabel('Source'); set(gcf,'color',[1 1 1]); colorbar; shading flat; axis tight
figure; plot(0:max_age,Lj); xlabel('Age'); ylabel('Length'); set(gcf,'color',[1 1 1]); box off; grid on
figure; plot(0:max_age,Wj); xlabel('Age'); ylabel('Weight'); set(gcf,'color',[1 1 1]);box off; grid on
end
Lj(1)=[]; Wj(1)=[];
Wij=repmat(Wj,numpatches,1); %weight at each age class, replicated across patches
Bij_initial=ones(numpatches,max_age).*(TargetBiomass/(numpatches*max_age));
Fsum=0; %no harvest
theta=0; %placeholder
% break
%% tune free parameters
if tune_parameters_1YES_0NO==1
    %initial guesses
    Rmax_guess=Rmax; 
    alphaCRguess=alpha1;
    for Ti=1:Tune_iters
        disp(['--Tuning iteration ',num2str(Ti),' --'])
        
  %     RMAX
        Rmax_LB=Rmax_guess/100;
        Rmax_UB=Rmax_guess*100;
        prob.obj = 'tuneRmax';
        prob.guess = Rmax_guess;
        prob.lo = Rmax_LB;
        prob.up = Rmax_UB;
        tic
    %     [FVAL,X, g, hess, inform, EXITFLAG]=npsol(prob); %optimize
        % [X,FVAL,EXITFLAG] = FMINCON(@tuneRmax,Rmax_guess,[],[],[],[],Rmax_LB,Rmax_UB)
        [X,FVAL,EXITFLAG] = fzero(@tuneRmax,Rmax_guess);
        disp(['Tuned Rmax = ',num2str(X),'; Took ',num2str(toc/60),' minutes. Flag = ',num2str(EXITFLAG),])
        Rmax_guess=X;
        Rmax=X;

        % AlphaCR
        alphaCR_LB=alpha1/100;
        alphaCR_UB=alpha1*100;
        prob.obj = 'tuneAlphaCR';
        prob.guess = alphaCRguess;
        prob.lo = alphaCR_LB;
        prob.up = alphaCR_UB;
        tic
    %     [FVAL,X, g, hess, inform, EXITFLAG]=npsol(prob); %optimize
        % [X,FVAL,EXITFLAG] = FMINCON(@tuneAlphaCR,alphaCRguess,[],[],[],[],alphaCR_LB,alphaCR_UB)
        [X,FVAL,EXITFLAG] = fzero(@tuneAlphaCR,alphaCRguess);
        disp(['Tuned AlphaCR = ',num2str(X),'; Took ',num2str(toc/60),' minutes. Flag = ',num2str(EXITFLAG),])
        alphaCRguess=X;
        alpha1=X;
        warning off MATLAB:divideByZero
        beta_i_initial=alpha1./(Rmax.*domain_suit_hab);
        warning on MATLAB:divideByZero
        
    end
    disp(['Final tuned values: Rmax = ',num2str(Rmax),'; alpha1 = ',num2str(alpha1),';'])
else %do not tune Rmax and alpha
    disp(['Manual values: Rmax = ',num2str(Rmax),'; alpha1 = ',num2str(alpha1)])
end
%Show tuned results
hide0_show1_figures=0;
MOPmodel2
disp(['CR = ',num2str(CR),])
disp(['CRgoal = ',num2str(CRgoal),])
KarryNumbers=sum(sum(Nij(:,age_mature:end)));
disp(['Carrying capacity in biomass        = ',num2str(K),])
disp(['Target Carrying capacity in biomass = ',num2str(TargetBiomass),])
disp(['Carrying capacity in numbers = ',num2str(KarryNumbers),])
Bij_initial=Bij; %set initial conditions
Kvirgin=K;
Bi_initial=sum(Bij_initial,2); %virgin biomass in each patch (for static model)

%THETA NOTES - INVERSE MODEL
% Revenue= pr*B*(1-1/(1+F+delta)*(F/(F+delta) %general formula
% MR=pr*B*(1-1/(1+F+delta)*(F/(F+delta) %as F goes to zero
% Cost = T*F
% MC = T*F as F goes to zero
% set MC=MR, cancel out F in the numerator:
% T=pr*B*(1-1/(1+F+delta)*(1/(F+delta)
% reduce B to theta_percent level, and set F to zero because interested in marginal rev/cost
% T=pr*B*Tpercent*(1-1/(1+0+delta)*(1/(0+delta))
% simplify
% T=pr*B*Tpercent*(1-1/(1+delta)*(1/(delta))
% theta=(pr*((Kvirgin/numpatches)*theta_percent))*(1-1/(1+delta))*(1/delta);

%THETA NOTES - EXP MODEL
% Revenue= pr*B*(1-exp(-(F+delta)))*(F/(F+delta) %general formula
% MR=pr*B*(1-exp(-(F+delta))*(F/(F+delta) %as F goes to zero
% Cost = T*F
% MC = T*F as F goes to zero
% set MC=MR, cancel out F in the numerator:
% T=pr*B*(1-exp(-(F+delta)))*(1/(F+delta)
% reduce B to theta_percent level, and set F to zero because interested in marginal rev and cost
% T=pr*B*Tpercent*(1-exp(-(0+delta)))*(1/(0+delta))
% simplify
% T=pr*B*Tpercent*(1-exp(-delta))*(1/(delta))
theta=pr*((Kvirgin/numpatches)*theta_percent)*(1-exp(-delta))*(1/(delta));
% break
%% Iterative search for optimal Fsum w/o NEW MPAs
MPA_locs(existMPAs)=1; %close NCI MPAs
NCI_MPA_locs=MPA_locs;
tic
hide0_show1_figures=0;
LB=0;
UB=Fsum_UB_prop_wrt_sumFishingHab*sum(domain_suit_hab(MPA_locs==0)); %.3*numpatches; 
Fsum_range=[LB:(UB/Fsum_NumIters):UB];
Fsum_i=0;
for Fsum=Fsum_range
    disp([,num2str(length(Fsum_range)-Fsum_i),])
    Fsum_i=Fsum_i+1;
    MOPmodel2
    Profit_y=sum(Profit_iy,1);
    Profit_fleet(Fsum_i)=Profit_y(end);
    Bijy_fleet(Fsum_i)=sum(sum(Bijy(:,:,end)));
    Yy=sum(Yiy,1);
    Yield_fleet(Fsum_i)=Yy(end);
end
disp(['Iterative search took ',num2str(toc/60),' minutes'])
loc=find(Profit_fleet==max(Profit_fleet));
disp(['Iterative search: Max profit is $',num2str(max(Profit_fleet)),' at Fsum = ',num2str(Fsum_range(loc)),'. Yield = ',num2str(Yield_fleet(loc)),])

figure
plot(Fsum_range,Yield_fleet,'k');
xlabel('$\bar{F}$','Interpreter','latex');
ylabel('Yield [kg]');
title('Fishery')

figure
plot(Fsum_range,Profit_fleet,'k');
xlabel('$\bar{F}$','Interpreter','latex');
ylabel('Profit [$]');
title('Fishery')
grid on

if stop_before_optimization==1
    break
    disp('stop_before_optimization==1')
end
% break
%% Find optimal Fsum w/o NEW MPAs
tic
disp('---Find optimal Fsum')
[I,J]=find(Profit_fleet==max(Profit_fleet));
hide0_show1_figures=0;
x0=Fsum_range(J); %Fsum_range(loc); %set guess at iterative max (determined above) 
% LB=0;
% UB=0.3*numpatches; %max(Fsum_range);
tic
%fmincon==================================================
% warning off
% options_OptFsum=optimset('Display','off');
% [Fsum_opt NEG_FleetProfit_Fsum_opt flag]=fmincon(@OptFsum,x0,[],[],[],[],LB,UB,[],options_OptFsum);
% warning on
%npsol==================================================
    TolX_value=1e-4*numpatches; %Set tolerance for Fsum search. 1e-4 is the default TolX. 
    TolFun_value=TolX_value;
options_fmincon_OptFsum=optimset('Algorithm','interior-point','Display','off','TolX',TolX_value,'TolFun',TolFun_value);
prob.obj = 'OptFsum';
prob.guess = [x0];
prob.lo = [LB];
prob.up = [Fsum_UB_prop_wrt_sumFishingHab*sum(domain_suit_hab(MPA_locs==0))];
% [NEG_FleetProfit_Fsum_opt,Fsum_opt, g, hess, flag, iter]=npsol(prob);  % OPTIMIZATION
[Fsum_opt,NEG_FleetProfit_Fsum_opt,flag,output]=fmincon(@OptFsum,prob.guess,[],[],[],[],prob.lo,prob.up,[],options_fmincon_OptFsum); %v7
%=========================================================

disp(['Optimization took ',num2str(toc/60),' minutes'])
FleetProfit_Fsum_opt=-NEG_FleetProfit_Fsum_opt;
%find biomass at optimum Fsum
Fsum=Fsum_opt;
hide0_show1_figures=0;
MOPmodel2
Profit_i_Fsumopt_NoMPAs=Profit_iy(:,end); %profit in each patch w/o MPAs (what you would expect from a survey) - for static model
Fi_Fsum_opt=Fi;
Bijy_Fsum_opt=sum(sum(Bijy(:,:,end)));
Bi_Fsum_opt=sum(Bijy(:,:,end),2); %biomass in each patches w/o MPAs (what you would expect from a survey) - for static model
disp(['Fsum_opt = ',num2str(Fsum_opt),' Fi = ',num2str(min(Fi_Fsum_opt)),'-',num2str(max(Fi_Fsum_opt)),'; FleetProfit_Fsum_opt = ',num2str(FleetProfit_Fsum_opt),'; Exit flag = ',num2str(flag),])

figure
plot(domain_suit_hab(existMPAs),Bi_Fsum_opt(existMPAs),'bo');
xlabel('Habitat [km^2]')
ylabel('Biomass [kg]')
% title('NCI MPAs, given optimal Fsum')
hold on

X_existMPAs=domain_suit_hab(existMPAs);
X_existMPAs_log=log(X_existMPAs);
Y_existMPAs=Bi_Fsum_opt(existMPAs);
Log_adder=1; %because some habitat values are less than 1
% polyfit_out=polyfit(X_existMPAs_log,Y_existMPAs,1); 
% [B_existMPAs_regress,BINT,R,RINT,STATS] = regress(Y_existMPAs,[(log(X_existMPAs)), ones(length(X_existMPAs),1)]); %with nonzero intercept
% [B_existMPAs_regress,BINT,R,RINT,STATS] = regress(Y_existMPAs,log(X_existMPAs)); B_existMPAs_regress(2)=0; %force intercept to be zero
[B_existMPAs_regress,BINT,R,RINT,STATS] = regress(Y_existMPAs,log(X_existMPAs+Log_adder)); B_existMPAs_regress(2)=0; %NEW force intercept to be zero
X_existMPAs_interp=[floor(min(X_existMPAs)):ceil(max(domain_suit_hab))];
%use this equation to interp biomass in other patches if they wereconverted to MPAs
% Y_existMPAs_interp=polyfit_out(1).*log(X_existMPAs_interp)+polyfit_out(2); 
% Y_existMPAs_interp=B_existMPAs_regress(1).*log(X_existMPAs_interp)+B_existMPAs_regress(2); 
Y_existMPAs_interp=B_existMPAs_regress(1).*log(X_existMPAs_interp+Log_adder)+B_existMPAs_regress(2); %NEW
plot(X_existMPAs_interp,Y_existMPAs_interp,'k-')

%static manager uses this to estimate biomass in patches if they were closed
% Kunfished_hat=polyfit_out(1).*log(domain_suit_hab)+polyfit_out(2); %based on polyfit
Kunfished_hat=B_existMPAs_regress(1).*log(domain_suit_hab+Log_adder)+B_existMPAs_regress(2); %based on regress
% Kunfished_hat(Kunfished_hat<0)=0; %no negative biomass levels. NEW: commented out because there shouldn't be any neg values in the first place
Kunfished_hat(existMPAs)=Bi_Fsum_opt(existMPAs); %manager know biomass in NCIs  so make exact estimates there
% note: the above strategy should work pretty well for the static manager
% when s/he closes just a few more patches. But, when many patches are
% closed they act synergistically to elevate biomass greater than what
% would be predicted based on the regression using just the 12 NCI MPAs 
tmp1=domain_suit_hab; tmp1(existMPAs)=[]; %just look at non-MPA patches
tmp2=Kunfished_hat; tmp2(existMPAs)=[]; %just look at non-MPA patches
plot(tmp1,tmp2,'k^');
legend('11 Existing MPAs','Regression','Interpolation','location','southeast')
legend('11 Existing MPAs',['Fit (R^2 = ' num2str(STATS(1)) ')'],'Interpolation','location','southeast')
title('Static manager estimate')
box off
set(gcf,'color',[1 1 1]);
plot(domain_suit_hab(existMPAs),Bi_Fsum_opt(existMPAs),'bo');

figure
plot(Bi_initial,Kunfished_hat,'bo')
xlabel('A world without fishing')
ylabel('Interpolation of particular closures')
title('Biomass in each patch')

% figure
% plot(domain_suit_hab(existMPAs),Bi_Fsum_opt(existMPAs)./domain_suit_hab(existMPAs),'bo');
% xlabel('Habitat [km^2]')
% ylabel('Density [kg/km^2]')
% title('NCI MPA patches, given optimal Fsum')
% break
%% Net present value
disp('-----------------------NPV')
%COMMENT OFF TO DO NPV (OTHERWISE EQUILIBRIUM RESULTS):
NPVon=0; 
if NPVon==1
    Bij_initial=Bij; %set initial conditions based on optimal no-new MPA management (just solved above)
    max_simulation_length=30;
    discount_rate=0.05;
    discounted_value_y=1./((1:max_simulation_length).^discount_rate);
    figure
    plot(1:max_simulation_length,discounted_value_y)
    xlabel('Future year')
    ylabel('Discounted value')
    set(gcf,'color',[1 1 1]);
    box off
    tmp1=NaN(numpatches, max_age, max_simulation_length);
    discounted_value_iy=repmat(discounted_value_y,numpatches,1);%for multiplying by Profit_iy
    for tmp2=1:length(tmp1(1,:,1));
        tmp1(:,tmp2,:)=discounted_value_iy;
    end
    discounted_value_ijy=tmp1; %for multiplying by Bijy

    %rerun no-new MPA results so can use in scaling below results
    Fsum=Fsum_opt;
    MOPmodel2
    tmp1=sum(Profit_iy,1); 
    FleetProfit_Fsum_opt=tmp1(end);
    %for the static model:
    Profit_i_Fsumopt_NoMPAs=Profit_iy(:,end); %profit in each patch w/o MPAs (what you would expect from a survey) - for static model
    % Fi_Fsum_opt=Fi;
    Bijy_Fsum_opt=sum(sum(Bijy(:,:,end)));
    Bi_Fsum_opt=sum(Bijy(:,:,end),2); %biomass in each patches w/o MPAs (what you would expect from a survey) - for static model
    disp(['NPV: Fsum_opt = ',num2str(Fsum_opt),' Fi = ',num2str(min(Fi_Fsum_opt)),'-',num2str(max(Fi_Fsum_opt)),'; FleetProfit_Fsum_opt = ',num2str(FleetProfit_Fsum_opt),'; Exit flag = ',num2str(flag),])

    figure
    plot(domain_suit_hab(existMPAs),Bi_Fsum_opt(existMPAs),'bo');
    xlabel('Habitat [km^2]')
    ylabel('Biomass [kg]')
    % title('NCI MPAs, given optimal Fsum')
    hold on

    X_existMPAs=domain_suit_hab(existMPAs);
    X_existMPAs_log=log(X_existMPAs);
    Y_existMPAs=Bi_Fsum_opt(existMPAs);
    Log_adder=1; %because some habitat values are less than 1
    % polyfit_out=polyfit(X_existMPAs_log,Y_existMPAs,1); 
    % [B_existMPAs_regress,BINT,R,RINT,STATS] = regress(Y_existMPAs,[(log(X_existMPAs)), ones(length(X_existMPAs),1)]); %with nonzero intercept
    % [B_existMPAs_regress,BINT,R,RINT,STATS] = regress(Y_existMPAs,log(X_existMPAs)); B_existMPAs_regress(2)=0; %force intercept to be zero
    [B_existMPAs_regress,BINT,R,RINT,STATS] = regress(Y_existMPAs,log(X_existMPAs+Log_adder)); B_existMPAs_regress(2)=0; %NEW force intercept to be zero
    X_existMPAs_interp=[floor(min(X_existMPAs)):ceil(max(domain_suit_hab))];
    %use this equation to interp biomass in other patches if they wereconverted to MPAs
    % Y_existMPAs_interp=polyfit_out(1).*log(X_existMPAs_interp)+polyfit_out(2); 
    % Y_existMPAs_interp=B_existMPAs_regress(1).*log(X_existMPAs_interp)+B_existMPAs_regress(2); 
    Y_existMPAs_interp=B_existMPAs_regress(1).*log(X_existMPAs_interp+Log_adder)+B_existMPAs_regress(2); %NEW
    plot(X_existMPAs_interp,Y_existMPAs_interp,'k-')

    %static manager uses this to estimate biomass in patches if they were closed
    % Kunfished_hat=polyfit_out(1).*log(domain_suit_hab)+polyfit_out(2); %based on polyfit
    Kunfished_hat=B_existMPAs_regress(1).*log(domain_suit_hab+Log_adder)+B_existMPAs_regress(2); %based on regress
    % Kunfished_hat(Kunfished_hat<0)=0; %no negative biomass levels. NEW: commented out because there shouldn't be any neg values in the first place
    Kunfished_hat(existMPAs)=Bi_Fsum_opt(existMPAs); %manager know biomass in NCIs  so make exact estimates there
    % note: the above strategy should work pretty well for the static manager
    % when s/he closes just a few more patches. But, when many patches are
    % closed they act synergistically to elevate biomass greater than what
    % would be predicted based on the regression using just the 12 NCI MPAs 
    tmp1=domain_suit_hab; tmp1(existMPAs)=[]; %just look at non-MPA patches
    tmp2=Kunfished_hat; tmp2(existMPAs)=[]; %just look at non-MPA patches
    plot(tmp1,tmp2,'k^');
    legend('11 Existing MPAs','Regression','Interpolation','location','southeast')
    legend('11 Existing MPAs',['Fit (R^2 = ' num2str(STATS(1)) ')'],'Interpolation','location','southeast')
    title('NPV: Static manager estimate')
    box off
    set(gcf,'color',[1 1 1]);
    plot(domain_suit_hab(existMPAs),Bi_Fsum_opt(existMPAs),'bo');

    figure
    plot(Bi_initial,Kunfished_hat,'bo')
    xlabel('A world without fishing')
    ylabel('Interpolation of particular closures')
    title('NPV: Biomass in each patch')

    Fsum=-1; %no harvest
    MOPmodel2
    tmp1=squeeze(sum(Bijy,1));
    tmp2=sum(tmp1,1);
    Kvirgin=tmp2(end);

end
% break
%% STATIC1: Optimally design MPAs wrt profit-biomass tradeoff - % MPA NOT SPECIFIED
STATIC1_v1; %FULLY STATIC *********PUT THIS BEFORE DYNAMIC MODEL, AND FEED RESULTS ************
% break
%% DYNAMIC: Optimally design MPAs wrt profit-biomass tradeoff - % MPA NOT SPECIFIED
%scaler for tracing EF (0=all conservation; 1=all payoff)
% Payoff_Conservation_EF_alpha_dynamic_range=[0:0.2:1]; %test runs
% Payoff_Conservation_EF_alpha_dynamic_range=[0:0.1:1]; %lo res
% Payoff_Conservation_EF_alpha_dynamic_range=[0:0.05:1];  %34 HOURS
% Payoff_Conservation_EF_alpha_dynamic_range=[0:0.01:1]; %Total hours of entire EF heuristic: 181.8204
Payoff_Conservation_EF_alpha_dynamic_range=Payoff_Conservation_EF_alpha_STATIC1_range;
% ^^^^ USE WHATEVER SIZE VECTOR CHOSEN BY STATIC MODEL

%shells
Payoff_Conservation_EF_wrt_alpha=NaN(2,length(Payoff_Conservation_EF_alpha_dynamic_range)); %storing payoff and conservation outcomes
tictoc_wrt_Payoff_Conservation_EF_wrt_alpha=NaN(1,length(Payoff_Conservation_EF_alpha_dynamic_range)); %time keeper
MPAlocs_Payoff_Conservation_EF_wrt_alpha=NaN(length(domain_suit_hab),length(Payoff_Conservation_EF_alpha_dynamic_range)); %for storing MPA locs solns
BestObjective_wrt_alpha=NaN(1,length(Payoff_Conservation_EF_alpha_dynamic_range));
PointHistory_wrt_alpha=NaN(1,length(Payoff_Conservation_EF_alpha_dynamic_range));
PointObjectiveHistory_wrt_alpha=NaN(1,length(Payoff_Conservation_EF_alpha_dynamic_range));
func_evals_wrt_alpha=NaN(1,length(Payoff_Conservation_EF_alpha_dynamic_range));
Fi_EF_wrt_alpha=NaN(length(domain_suit_hab),length(Payoff_Conservation_EF_alpha_dynamic_range));
EXITFLAG_dynamic=NaN(1,length(Payoff_Conservation_EF_alpha_dynamic_range));

%fill shell with extreme (alpha=0) (note: don't know result for alpha=1 yet!)
%note: assuming that Fsum opt produces optimal MPA design for maximizing profit (may be 0% MPAs)
Payoff_Conservation_EF_wrt_alpha(1,1)=0; %payoff with alpha=0
Payoff_Conservation_EF_wrt_alpha(2,1)=Kvirgin; %conservation with alpha=0
% Payoff_Conservation_EF_wrt_alpha(1,end)=FleetProfit_Fsum_opt; %payoff with alpha=1
% Payoff_Conservation_EF_wrt_alpha(2,end)=Bijy_Fsum_opt; %conservation with alpha=1
MPAlocs_Payoff_Conservation_EF_wrt_alpha(:,1)=1; %MPA design with alpha=0
% MPAlocs_Payoff_Conservation_EF_wrt_alpha(:,end)=(Fi_Fsum_opt==0); %MPA design with alpha=1
Fi_EF_wrt_alpha(:,1)=0; %no harvest bc all MPAs


% %fill shell with extreme (alpha=0, thus 100% MPA)
% %note: for alpha=1 we do not know if MPAs are involved or not (or, could
% %base on Fsum opt? <- But MPA plan may be equivalent...)
% Payoff_Conservation_EF_wrt_alpha(1,1)=FleetProfit_Fsum_opt; %payoff with 0% MPA
% Payoff_Conservation_EF_wrt_alpha(2,1)=Bijy_Fsum_opt; %conservation with 0% MPA
% Payoff_Conservation_EF_wrt_alpha(1,end)=0; %payoff with 100% MPA
% Payoff_Conservation_EF_wrt_alpha(2,end)=Kvirgin; %conservation with 100% MPA
% MPAlocs_Payoff_Conservation_EF_wrt_alpha(:,1)=0; %0=no MPA in that loc
% MPAlocs_Payoff_Conservation_EF_wrt_alpha(:,end)=1; %1=yes MPA in that loc

%now find payoff-conservation EF for every 0<alpha<1  
x0_Fsumguess=Fsum_opt; %use as starter, but then update wrt alpha and MPAlocs
Payoff_Conservation_EF_alpha_dynamic_i=1;
tic
index=0;
for Payoff_Conservation_EF_alpha_dynamic=Payoff_Conservation_EF_alpha_dynamic_range(2:end);
% for Payoff_Conservation_EF_alpha_dynamic=Payoff_Conservation_EF_alpha_dynamic_range(17:end);    
    Payoff_Conservation_EF_alpha_dynamic_i=Payoff_Conservation_EF_alpha_dynamic_i+1;
    index=index+1;
    
    disp(['--Payoff_Conservation_EF_alpha_dynamic = ',num2str(Payoff_Conservation_EF_alpha_dynamic),])
%     Fsum=Fsum_opt; %that is if do not optimize Fsum each time
    
    Dims=numpatches-length(existMPAs); % #patches not including those fixed as MPAs
    if index==1
        MaxSingleChangeTrials=1e4; %default 1e4
        TabuOnOff=1; %default=1
        %Crow: have it start at reasonable #MPAs:
        Point = zeros(Dims,1); % Initiate point
        idx = randperm(Dims); idx = idx( 1:(floor((1-Payoff_Conservation_EF_alpha_dynamic)*Dims)));
        Point( idx ) = 1; Point_initial=Point;
    end %otherwise update from previous iteration
    
    t1start = tic;
%     [ BestPoint, BestObjective, PointHistory, PointObjectiveHistory, func_evals ]=...
%         BinaryTabuMonteCarlo(@Payoff_Conservation_EFvalue_v1, Dims, MaxSingleChangeTrials,TabuOnOff, Point_initial);
    
    % ga
   options_ga=gaoptimset('Display','off','InitialPopulation', Point_initial');
    [X,FVAL,EXITFLAG,OUTPUT]=ga(@Payoff_Conservation_EFvalue_v1,Dims,[],[],[],[],zeros(1,Dims),ones(1,Dims),[],[1:Dims],options_ga);
     BestPoint=X';
     BestObjective=FVAL;
    func_evals=OUTPUT.funccount;
    %%%
    
    Point_initial=BestPoint; %update heuristic starting point
    tictoc_wrt_Payoff_Conservation_EF_wrt_alpha(Payoff_Conservation_EF_alpha_dynamic_i)=toc(t1start)/60;
    BestObjective_wrt_alpha(Payoff_Conservation_EF_alpha_dynamic_i)=BestObjective;
%     PointHistory_wrt_alpha(Payoff_Conservation_EF_alpha_dynamic_i)=PointHistory;
%     PointObjectiveHistory_wrt_alpha(Payoff_Conservation_EF_alpha_dynamic_i)=PointObjectiveHistory;
    func_evals_wrt_alpha(Payoff_Conservation_EF_alpha_dynamic_i)=func_evals;
    EXITFLAG_dynamic(Payoff_Conservation_EF_alpha_dynamic_i)=EXITFLAG;

    tmp=NCI_MPA_locs'; %set vector just with NCI MPAs
    tmp(NCI_MPA_locs==0)=BestPoint; %fill itn the new MPAs by heuristic
    BestPoin_allMPA_locs=tmp;
    MPAlocs_Payoff_Conservation_EF_wrt_alpha(:,Payoff_Conservation_EF_alpha_dynamic_i)=BestPoin_allMPA_locs; %store MPA locs soln
    
    %run output MPA design and store payoff and conservation outputs
    MPA_locs=BestPoin_allMPA_locs'; %locs that are fished (0) or MPAs (1) -  (input from Heuristic)
    MOPmodel2
    tmp1=sum(Profit_iy,1);
    tmp2=tmp1(end); 
    Payoff_Conservation_EF_wrt_alpha(1,Payoff_Conservation_EF_alpha_dynamic_i)=tmp2; %payoff
    tmp1=sum(sum(Bijy(:,:,end)));
    Payoff_Conservation_EF_wrt_alpha(2,Payoff_Conservation_EF_alpha_dynamic_i)=tmp1; %conservation
    Fi_EF_wrt_alpha(:,Payoff_Conservation_EF_alpha_dynamic_i)=Fi;
    data_save_midrun=Payoff_Conservation_EF_wrt_alpha;
    save data_save_midrun
    disp(['#MPAs=',num2str(length(find(MPA_locs==1))),'; Fsum=',num2str(Fsum),'; x0_Fsumguess-Fsum_out=',num2str(x0_Fsumguess-Fsum),'; took ',num2str(tictoc_wrt_Payoff_Conservation_EF_wrt_alpha(Payoff_Conservation_EF_alpha_dynamic_i)),' minutes; EXITFLAG=',num2str(EXITFLAG)'])
    disp(['Payoff_Conservation_EF_alpha_dynamic=',num2str(Payoff_Conservation_EF_alpha_dynamic),'; out (goal is NEG ONE)=',num2str(BestObjective),'. Saved as data_save_midrun' ]);
end
disp(['Minutes of each heuristic: ',num2str(tictoc_wrt_Payoff_Conservation_EF_wrt_alpha),])
disp(['Total hours of entire EF heuristic: ',num2str(nansum(tictoc_wrt_Payoff_Conservation_EF_wrt_alpha)/60),])

figure %absolute outputs
plot(Payoff_Conservation_EF_wrt_alpha(2,:),Payoff_Conservation_EF_wrt_alpha(1,:),'b.')
xlabel('Biomass [kg]')
ylabel('Profit [$]')
set(gcf,'color',[1 1 1]); box off
title('Optimal MPA-Fleet outcomes')

%scale outputs
Payoff_Conservation_EF_wrt_alpha_scaled=Payoff_Conservation_EF_wrt_alpha;
Payoff_Conservation_EF_wrt_alpha_scaled(1,:)=Payoff_Conservation_EF_wrt_alpha_scaled(1,:)./FleetProfit_Fsum_opt; %scale payoff to fleet model max (note: not necessarily 0% MPA)
Payoff_Conservation_EF_wrt_alpha_scaled(2,:)=Payoff_Conservation_EF_wrt_alpha_scaled(2,:)./Kvirgin; %scale conseration to 100% MPA outputs

%indicate alpha next to point
figure %scaled outputs
plot(Payoff_Conservation_EF_wrt_alpha_scaled(2,:).*100,Payoff_Conservation_EF_wrt_alpha_scaled(1,:).*100,'bo','markerfacecolor','b')
xlabel('Biomass [%]')
ylabel('Profit [%]')
set(gcf,'color',[1 1 1]); box off
title('Optimal MPA-Fleet outcomes')
ylim([0 110])
for tmp=1:length(Payoff_Conservation_EF_alpha_dynamic_range)
    text_offset_x=.5;
    text_offset_y=1;
    text(Payoff_Conservation_EF_wrt_alpha_scaled(2,tmp).*100+text_offset_x,Payoff_Conservation_EF_wrt_alpha_scaled(1,tmp).*100+text_offset_y,num2str(Payoff_Conservation_EF_alpha_dynamic_range(tmp)))
end

data_save=Payoff_Conservation_EF_wrt_alpha_scaled;
save data_save

% break
% %% STATIC1: Optimally design MPAs wrt profit-biomass tradeoff - % MPA NOT SPECIFIED
% STATIC1_v1; %FULLY STATIC *********PUT THIS BEFORE DYNAMIC MODEL, AND FEED RESULTS ************
%% compare results
figure %scaled outputs
plot(Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(2,:).*100,Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(1,:).*100,'r^','markerfacecolor','r')
hold on
plot(Payoff_Conservation_EF_wrt_alpha_scaled(2,:).*100,Payoff_Conservation_EF_wrt_alpha_scaled(1,:).*100,'bs','markerfacecolor','b')
xlabel('Biomass [%]')
ylabel('Profit [%]')
set(gcf,'color',[1 1 1]); box off
legend('Informed static model','Dynamic model')
ylim([0 110])

% x=Payoff_Conservation_EF_wrt_alpha_scaled(2,:).*100;
% y=Payoff_Conservation_EF_wrt_alpha_scaled(1,:).*100;
% dt=DelaunayTri(x',y');
% k=convexHull(dt);
% plot(x(k), y(k), 'r')

%compare static and dynamic MPA designs
%any zeros indicate not change
%pos=patch was dynamic MPA but static fished
%neg=patch was dynamic fished by static MPA
MPAdiff=MPAlocs_Payoff_Conservation_EF_wrt_alpha-MPAlocs_Payoff_Conservation_EF_wrt_alpha_STATIC1;
figure
surface(Payoff_Conservation_EF_alpha_dynamic_range,1:numpatches,MPAdiff)
shading flat
axis tight
xlabel('\alpha')
ylabel('Patch')
set(gcf,'color',[1 1 1]);
% colorbar
% title('1=dynamic MPA but static fished; -1=dynamic fished but static MPA')
title('Red=dynamic MPA only; Blue= static MPA only; Green=no difference')

tmp1=MPAdiff; 
tmp2=tmp1>0; tmp3=sum(tmp2,1);
tmp4=tmp1<0; tmp5=sum(tmp4,1);
tmp6=tmp3+tmp5;
figure
plot(Payoff_Conservation_EF_alpha_dynamic_range,tmp3, 'r-')
hold on
plot(Payoff_Conservation_EF_alpha_dynamic_range,tmp5, 'b-')
plot(Payoff_Conservation_EF_alpha_dynamic_range,tmp6, 'k-')
xlabel('\alpha')
ylabel('No. patches')
legend('Dynamic MPA but static fished','Dynamic fished but static MPA','Mismatch either way')
axis tight
set(gcf,'color',[1 1 1]);
box off

% #MPAs per alpha
figure
plot(Payoff_Conservation_EF_alpha_dynamic_range,sum(MPAlocs_Payoff_Conservation_EF_wrt_alpha_STATIC1,1),'r^')
hold on
plot(Payoff_Conservation_EF_alpha_dynamic_range,sum(MPAlocs_Payoff_Conservation_EF_wrt_alpha,1),'bs')
xlabel('\alpha')
ylabel('No. MPAs')
legend('Static','Dynamic')
axis tight
set(gcf,'color',[1 1 1]);
box off
%%
if NPVon==0
    data_save=Payoff_Conservation_EF_wrt_alpha_scaled;
    save data_save_equil
elseif NPVon==1
    data_save_NPV=Payoff_Conservation_EF_wrt_alpha_scaled;
    save data_save_NPV
end
    
break
VoDI_figs_v1

